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Abstract 

In spite of many attempts to model dense granular flow, there is still no general theory capable of 
describing different types of flows, such as gravity- driven drainage in silos and wall-driven shear flows 
in Couette cells. Here, we summarize our recent proposal of the Stochastic Flow Rule (SFR), which is 
able to describe these cases in good agreement with experiments, and we focus on testing the theory 
in more detail against brute- force simulations with the discrete-element method (DEM). The SFR is a 
general rate- independent constitutive law for plastic flow, based on diffusing "spots" of fluidization. In 
the case of quasi-2D granular materials, we assume limit-state stresses from Mohr-Coulomb plasticity and 
postulate that spots undergo biased random walks along slip-lines, driven by local stress imbalances. We 
compare analytical predictions of the SFR against DEM simulations for silos and Couette cells, carrying 
out several parametric studies in the latter case, and find good agreement. 

1 Introduction 

Even after centuries of study, dense granular flow remains an enigma for both physicists and engineers [121 E] • 
For dilute, collisional granular media, the kinetic theory of dissipative gases has been quite successful at 
describing the mean flow statistics \J_^ . On the opposite side of the spectrum, where the material is solid-like 
and possesses a yield function, the stresses and (to a lesser extent) the dynamics can be adequately described 
by methods of modern soil mechanics [2F^'33j. However, the intermediate regime where the flow is slow, 
steady, and dense, has resisted many theoretical attempts PfM]. 

We have recently proposed the Stochastic Flow Rule (SFR) for dense granular flows [14J. It is a rate- 
independent constitutive law for steady flow which we believe naturally extends the applicability of traditional 
soils plasticity. Since it is rooted in a general theory for granular stresses, the SFR can be applied in any 
environment with two-dimensional (2D) symmetry, a level of generality not common in other models. Though 
the mean stress field is assumed to obey a continuum law, the SFR does not, however, treat the flow as 
continuous, but rather as the superposition of many discrete plastic flow events. The carriers of plastic flow 
are "spots" of local fluidization, which diffuse through the material in response to stress imbalances. The 
Spot Model provides a precise and general mechanism for random-packing dynamics [4J , which can be either 
used as a basis for multiscale particle simulations [30j or averaged to obtain simple continuum flow equations, 
which can be solved exactly. Here we take the latter approach and compare the analytical predictions of the 
SFR to brute- force simulations of dense granular flows by the Discrete Element Method (DEM). 

We begin by summarizing the theory, starting with the notion of the "spot", as a discrete carrier of 
granular motion. We then review limit-state Mohr-Coulomb plasticity for stresses, and proceed to describe 
how the local stress state biases the motion of spots. With the SFR fully defined, we then solve for the mean 
flow field in the draining-silo and Couette annular-flow geometries and compare to DEM simulations fo both 




Figure 1: Spots as carriers of plastic deformation in granular materials. Cartoon of basic spot motion. A spot 
of local fluidization, carrying some free volume, moves to the upper right causing a cooperative displacement 
of particles, to the lower left, opposing the spot displacement, (b) In silo drainage, spots are injected at the 
orifice and perform random walks biased upward by gravity, causing downward motion of particles. 



geometries. For annular fiow in particular, we test predictions of the dependence of the fiow profile on the 
particle contact friction, rotation rate of the inner wall, and inner wall radius, which constitute rigorous tests 
of the theory. 

2 Theory 

2.1 Spots and spot dynamics 

At the outset, we must first discuss a particular property of granular fiows underscoring their discrete nature. 
Unlike a true continuum, a steady flow field for a granular media is actually only steady in a time-averaged 
sense. There are constant microscopic fluctuations occurring in the flow caused chiefly by geometric packing 
constraints: grains cannot pass through other grains, so frequently the particle motion must fluctuate from 
the mean to achieve flow. 

In experiments and simulations of silo drainage, a pronounced correlation length for these fluctuations 
was observed [8l[9l[30j. This length, typically 3d to 5(i, where d is the particle diameter, implies that a grain 
does not move independent of its surroundings. Rather, it travels in a collective fashion, as if dragging a set 
of nearest neighbors with it as it moves. 

Based on this evidence for correlated motion, Bazant constructed the Spot Model for random packing 
dynamics and applied it to the case of granular drainage [4J . The model is based on a collective microscopic 
mechanism for flow, where particles move to oppose the random motion of "spots", as depicted in Figure 
[H When a spot displaces, the nearby affected particles respond by displacing in a block-like fashion in the 
opposite direction. The spot can thus be viewed as effectively carrying free volume and would typically be 
associated with a slight increase in the local interstitial volume of the packing. As a flrst approximation, 
identical spots undergo independent random walks through the material, although interactions between spots 
(e.g. carried by the stress fleld) and statistical variations in spot size and shape could be considered as well. 

Each step of its random trajectory, a spot chooses the displacement Axg following a biased random 
process. The mean displacement, which is effectively the deterministic part of the motion, is given by 

(Ax,) = AtBi (1) 

for step duration At and drift vector Di. The random part of the motion can be described at the simplest 
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level by 

(Ax, (g) Ax,) = 2AtD2 (2) 

where D2 is the symmetric tensor of diffusion coefficients (co- variance matrix) measuring the random fluctu- 
ations in each direction. These coefficients suffice to construct a continuous Fokker-Planck (drift-diffusion) 
equation to approximate the spot concentration, 

^ + V • (Dip,) = (V V) : (D2P,) + . . . (3) 

where omitted correction terms (in the Kramers-Moyall expansion) depend on higher moments of the spot 
displacement [29j. For example, in a simple model of silo drainage (ressembling the classical kinematic [27 
and void [T71 [21] models), the drift vector may be chosen to point uniformly upward, since particles drift 
downward due to gravity, with a constant diffusivity tensor to obtain spot trajectories as sketched in Figure 

In steady-state flow, the spot concentration satisfies the time-independent drift-diffusion equation, 

V • (Dip,) = (V ® V) : (Dsp,) (4) 

where A : B = AijBij. Once solved, the mean particle velocity field u can be found by superposing the 
effects of all spots on all the particles. This is easiest to see as a convolution of some influence function 
K;(r, r'), representing the amount of influence a spot centered at has on a particle at r, and the negative 
spot flux vector 

u* = -Dip, + V-(D2p,). (5) 

Thus the particle velocity u is 

u = J w{r,r')u%r')dr\ (6) 

The influence function, which in simple language describes the spot's shape, should have a characteristic 
width of three to five particle diameters, so that spots match known correlation length data. In many 
situations, u* is close to u since the convolution with w tends only to smooth out sharp changes in the spot 
flux density. 

There are several ways to use the multiscale description of the Spot Model, as have been demonstrated 
in the case of silo drainage. Given the drift field Di, diffusion tensor field D2, and for extra clarity the spot 
influence function equations (|4|) and ([6]) could be solved analytically or numerically to give the mean 
flow field; the statistical trajectories of particles could also be constructed, thus linking particle diffusion to 
the diffusion of free volume [4]. The microscopic basis for the Spot Model, however, also allows for explicit 
multiscale simulations of all particles in a granular flow, alternating between mesoscopic spot displacements 
and cooperative particle motion; this approach yields remarkably realistic flowing packings (compared to 
DEM simulations) in silo drainage [30j . This suggests that the spot mechanism provides a robust description 
of how granular flow occurs and could apply in a more universal sense, beyond just the gravity driven 
drainage from a silo. 

The difficult question, if we want to extend the spot model to an arbitrary environment, is how can we 
derive the spot dynamics (or at least Di and D2) from mechanical principles? The intuition behind the 
simple spot dynamics in silo drainage has no clear extension to other situations, such as an annular Couette 
cell, where the flow is unaffected by gravity. Instead, the spot dynamics must be related to stresses in the 
material, not only due to body forces, but also due to forces transmitted from moving boundary tractions. 



2.2 Quasi-static granular stresses 

In a dense granular material, forces are transmitted across discrete contact networks, which can be rather 
complicated and varying at the same length scale as the velocity profile or even the container geometry. 
Even the simplest question of how forces are transmitted in static granular heaps has been the subject of 
considerable controversy [H [TOl [22l [35l [15] . Nevertheless, it seems reasonable that a continuum theory could 
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at least describe the mean stress state in a steady flow as a starting point for applying the SFR. Some of the 
random fluctuations in the force network could also be manifested in the random motion of spots coupled 
to mean stresses. 

Most common continuum theories for solid-like material that can plastically flow utilize some type of 
elasticity, generating the stress tensor T from a small elastic strain via T = ^[E^] where the fourth rank 
tensor function ^ is determined by the chosen elasticity formulation. Elliptic systems of equations, such as 
the laws of linear isotropic elasticity, have solutions for any given set of boundary tractions or displacements, 
and thus complete boundary data must be supplied. 

This is problematic for granular materials since the boundary stresses in a static assembly are frequently 
the quantity of interest. For example, the pressure distribution under a sand-pile or the wall stresses in a bin 
have both been studied extensively, as noted above. Static granular boundary conditions typically develop 
internally as a response to the construction proceedure and often cannot be fluctuated without plastically 
deforming the assembly. 

Hyperbolic laws of stresses have the beneflt of constraint: if the tractions along one boundary region 
are chosen, their effect can propagate through the bulk to other boundaries and place constraints on the 
allowable tractions there. Some non-linear elasticity laws have been shown to yield unique solutions under 
restrained boundary information [T3] , but it appears a different continuum approach may be more naturally 
suited for granular materials. 

The traditional approach from soil mechanics is not a strain-based stress formulation, but rather an 
analysis based more closely on a yield function. The Coulomb yield function deflnes failure whenever \r/a\ = 
ji = tan^ where r is the shear stress on a plane, a the normal stress, and ji is an internal friction coefficient 
unique to the material. This yield function akin to a cohesionless dry friction law seems the natural choice for 
a continuum description of granular materials — applying higher compression locks the grains more securely, 
and thus makes it harder to shear them. Our experience with granular piles also direct us to this yield 
function. If we pour sand from a narrow spout onto a rough table, grains build up on the surface of the 
sand cone until the cone angle reaches a certain value at which point grains always slide down the surface, 
analogous to how a frictional block resting on a board will start to slide if the board is tilted above some 
critical angle. 

In the Critical State Theory of soils ^33j , the stresses throughout a moving soil eventually converge onto a 
"critical state line" where pressure p and stress deviator q are linearly related. This locus strongly resembles 
the Drucker-Prager yield criterion, a simplifled Coulomb yield law for 3D. The deformation rate in such slow 
flows does not signiflcantly contribute to the stress state. 

For this and other reasons we too will adopt a rate-independent formulation for granular stresses, which 
should owe well to the slow flowing regime we are primarily concerned with. To connect the flowing to 
the static regime, a common assertion is that the yield criterion is being met at all points in the material 
whether or not the material is flowing yet [26j . When static material is brought to steady motion under this 
assertion, the locus of static stresses already resembles the critical state line and thus no new stress laws are 
needed. 

Materials obeying this traditional hypothesis are said to be at limit-state. Let us focus entirely on 
2D geometries (plane strain) complete with a 2D stress tensor T deflned in the plane of the flow. By 
requiring the yield criterion to be met at all points, the limit-state assumption implies the following constraint 
^ |To| = sin(/) tr T) where Tq is the deviatoric stress tensor and |A| = a/A : A. A simple way to uphold 
this constraint is to rewrite the stress fleld in terms of two stress parameters (the Sokolovskii variables [36]): 
p the average pressure, and the so-called "stress-angle" denoting the angle from the horizontal to the 
major principal stress. The components of the 2D stress tensor are then 

Txx = — p(l + sin(/)cos 2?/^) (7) 
Tyy = — p(l — sin^cos 2?/^), (8) 
Txy = Tyx = -p sin sin 2?/^ (9) 

from which it can be seen that p = — trT/2 and tan 27/; = 2Ti2/(Tii — T22). The convection-free 2D 
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momentum balance law V • T + Fbody = then reduces to the two variable system of hyperbolic PDEs 



(1 + sin (j) cos 2tlj)px — 2p sin <p sin 2ip ip^ + sin cj) sin 2?/^ py + 2p sin (j) cos 2ip t/j, 
sin (/) sin 2?/^ px + 2p sin (j) cos 2?/^ ^/^^^ + (1 — sin (j) cos 2tlj)py + 2p sin </) sin 2?/^ ?/^, 



'y 



'y 



= Fi 



body* 



body 



(10) 

(11) 



A simple Mohr's circle analysis shows that the directions along which the yield criterion is met, the slip-lines, 
form at the angles ip ^ e from the horizontal where e = 7r/4 — 0/2. 

2.3 The stochastic flow rule 

To obtain flow from these stress equations, we will describe an unorthodox, yet effective, "stochastic flow 
rule" (SFR) based on postulating spot dynamics driven by local stress imbalances [M]. Unlike classical 
flow rules in plasticity which assume an ideal continuous medium, the SFR reflects the discreteness and 
randomness inherent to a granular packing. In its simplest embodiment, the SFR provides a unique way of 
defining the spot drift Di and diffusivity D2 in an arbitrary plane strain environment. 

The SFR is based physically on the notion that spots move by sliding on average along the (Mohr- 
Coulomb) slip-lines. However, there are two crucial sources of randomness: The slip-lines themselves are 
blurred by statistical fluctuations [28j, but, more fundamentally, the spot must randomly decide between 
one of two slip-lines through each material point to slide in each step. Given the discrete nature of the 
packing at the spot scale (several particle diameters), it would be geometrically impossible for the spot to 
cause particle displacements along both slip-lines at once. 

In principle, the SFR is very general and could apply to any limit-state stress field (not only for a 
granular material), but to arrive at the simplest possible model for granular flow, with only one parameter, 
we make several bold assumptions. The first is that of isotropic spot diffusion D2 = D2I which may be 
justified to some extent by slip-line blurring. To determine the magnitudes of the drift and diffusion, we 
reason qualitatively. A spot can be viewed as a local region of internal plastic failure, which perturbs the 
contact stresses on neigboring material elements. Such a perturbation can incite adjacent material to fail 
and when it does, the spot effectively propagates to a new location. For a typical spot size Lg, the distance 
of propagation should thus be roughly Lg as well. Accordingly, we constrain all lengths referenced in the 
definitions of the drift and diffusion, equations ([T]) and ([2|), to uphold this notion, i.e. |Di| = L^/At, and 
D2 = Ll/2/Si. All that remains is to determine the drift direction dg = Di/ |Di|, which is where stresses 
enter the model. 

We derive the spot drift from local stress imbalances upon fluidization as follows (see [14J for more 
details). Define a material cell as the roughly diamond-shaped region encompassed by two intersecting pairs 
of slip-lines, separated by Lg. If a spot enters a cell, the value of the internal friction along the cell boundaries 
should drop from /i to the kinetic value fik causing a change in the shear tractions. This fluctuation may 
generate a net force on the cell, which takes the simple form. 



This force pushes the material one way, and consequently a spot should move the other. Since there is a 
slip-line field, however, a spot cannot move in any arbitrary direction. Instead, the motion opposing the net 
force must be projected onto slip-lines to obtain allowable motions. By projecting — Fnet onto the slip-lines 
and averaging, we find that the spot drift direction is given by 



With these assumptions, the SFR is specified enough for multiscale simulations, alternating between 
continuum stress computation, mesoscale spot random walks, and cooperative particle displacements. Al- 
ternatively, as a continuum approximation, the mean flow profile can be constructed by a two-step process: 



Fnet = (1 - M/c/m) (Fbody " COS^ Vp). 



(12) 





(13) 



(14) 
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1. Calculate the steady spot concentration the drift- diffusion equation, which now has the simplified form 

V • {d,ps) = yV^p,. (15) 

2. Compute the mean particle velocity by convolving the spot infiuence function with the reversed spot 
fiux, 

" = "ii / "'("'"'^ (d.(r')p.(r') - YVp,(r')) dr'. (16) 

As can be seen from above, the spot step duration At has no effect on the fiow profile, aside from a 
simple rescaling of the velocity. Together with the fact that equation (p!5|) is homogeneous in ps , we observe 
a consequence of the SFR common to other rate-independent models: Any velocity field which solves the 
SFR equations is unique only up to a multiplicative constant. 



3 Solving for Flow in Simple Geometries 

The SFR, with some theoretical exceptions as detailed in , may apply as a general law in 2D plasticity. 
We will focus on two canonical cases which represent two very different types of granular fiow: gravity- 
driven drainage from a silo and forced shear fiow in an annular Couette cell. We are not aware of any other 
model, continuum or discrete, which can describe both of these cases, even qualitatively, so this will be a 
stringent test of the SFR. We emphasize there will be no fitting parameters used. The parameters and 
Ls are independently measured material properties, and their values are assumed not to depend on the fiow 
environment. 



3.1 Silo 

For a wide, 2D silo with smooth walls and a fiat bottom, the stress balance equations can be solved analyt- 
ically, giving 

/ TT fg{Zm-z) 

2 1 + sm 

where z is the distance from the silo bottom, Zm is the distance from the bottom to the free surface, and fg 
is the material's weight density. Applying equation (p!2]) gives Fnet pointing uniformly downward and thus 
dg = z. To visualize the stress field (via its slip-lines) simultaneously with the spot drift vector, see Figure 
[2](a). As such, the spot concentration solves 

For the narrowest possible silo opening, we have the boundary condition that ps(x,0) = S{x). In line with 
our intuition on jammed states, this implies the opening must be at least one spot width for any fiow to 
occur. Solving with a Fourier transform gives [14J 

_ exp{-x^/4bz) ^^^^ 
V^Trbz 

for 6 = L5/2 and consequently, by substituting ps into equation ([5]), we have u* • z (x — ^^^^^=^^. 

Since the velocity gradient should not rapidly fiuctuate in the bulk fiow region we are interested in, 
we utilize u u*. Diffusive spreading of the downward velocity, as predicted here, is a well-documented 
phenomenon in silo drainage [27l [32l [HI [37l [20j . Utilizing the Ls dependence, the typical range of b values 
is predicted to be l.bd to 2.5d which compares quite well to the results of documented experiments on silo 
diffusion which, when accumulated, give a total b range of 1.3d to 3.5d. To our knowledge, there is no other 
theory to predict 6, even qualitatively, from basic mechanical principles, so this should be viewed as a first 
success of the SFR. 
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Figure 2: Slip-lines fields (from MCP) and the spot drift field (from the SFR) displayed for (a) a wide silo 
draining under gravity, and (b) shearing in an annular Couette cell (no gravity). 



3.2 Annular Couette cell 

A much stronger test of the SFR comes from its ability to adapt to a completely different flow environment, 
without changing any parameters. In annular Couette flow, material fills the region between two rough 
cylinders and is sheared by rotating the inner cylinder while holding the outer stationary. To solve for the 
MCP stresses, we first convert the stress balance equations to polar coordinates (r, 0) and require that p and 
obey radial symmetry. This simplifies to 

9?/;* sin 2?/;* df] 2sin0 

dr r(cos 2?/;* + sine/)) ' dr r(cos 2?/^* + sin (/)) 

where rj = \ogp and ip* = ip -\- ^ — 0. We solve these equations numerically using fully rough inner wall 
boundary conditions tp*{rw) = f ~ ^ and any arbitrary value for r]{rw)' We find that p decreases with 
distance from the inner wall, implying Fnet points radially outward. Applying equations (p!3|) and (p!4j) gives 
a drift field that points inward, but gradually opposes the motion of the inner wall. 

Below, we solve for u* directly by enforcing that the flow point along 0. Since the material near the wall 
shears rapidly, we cannot assert u u* and must solve the full convolution, equation (|6]). For simplicity, 
we assume circular spots w{r^r') ex H[kf- — |r — r'|) where is a step function. (Velocity correlations 
measurements in silo experiments suggest that the spot shape may be anisotropic due to gravity, or more 
generally, local stresses, but this is not a major effect [14 .) To properly describe particle motion near the 
wall, we must decide how spots interact with walls so that the convolution integral can be carried out all 
the way to the inner wall. A simple claim is that u* must equal the wall velocity wherever spots overlap the 
wall — a no-slip condition of sorts. After convolution, the net effect of this assertion is that u is slightly flatter 
than u* close to the inner wall. We shall refer to this assumption about wall interactions as "wall hypothesis 
1" . We note that this is slightly different than the hypothesis used in [14] which claims the region within the 
walls can be viewed as a bath of non-diffusive spots which cause particles to move in a manner that mimics 
the rigid wall motion. We shall refer to this as "wall hypothesis 2" . Though our theory is flexible in terms 
of near wall specificities, the bulk behavior, regardless of which hypothesis is selected, is always the same. 
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4 Comparing SFR Predictions to Discrete Element Simulations 



4.1 Simulation Methods 

To test the above models, we carried out DEM simulations using a modified version of the Large-scale 
Atomic/ Molecular Massively Parallel Simulator (LAMMPS) developed at Sandia National Laboratories 
Following a similar approach to previous work by our group [30', "ST] and others \3^, T5] we considered 
simulations of spherical particles interacting using a modified version of the model developed by Cundall and 
Strack [TT] to model cohesionless particulates. Monodisperse spheres with diameter d interact according to 
Hertzian, history dependent contact forces. For a distance r between a particle and its neighbor, when the 
particles are in compression, so that (5 = d — |r| > 0, then the two particles experience a force F = + F^, 
where the normal and tangential components are given by 

Fn = [knSn - ^) , F, = [-hAst - ^) . (21) 

Here, n = r/ |r|. and are the normal and tangential components of the relative surface velocity, and kn^t 
and 7n,t are the elastic and viscoelastic constants, respectively. As^ is the elastic tangential displacement 
between spheres, obtained by integrating tangential relative velocities during elastic deformation for the 
lifetime of the contact, and is truncated as necessary to satisfy a local Coulomb yield criterion |Ft| < 
/ic |Fn|. Particle- wall interactions are treated identically, though the particle- wall friction coefficient /j.^ is 
set independently. 

For the current simulations we set kt = jkn, and choose kn = '^x lO^mg/d. This is significantly less than 
would be realistic for many hard granular materials, such as glass beads, where we expect kn > lO^^mg/d^ 
but such a spring constant would be prohibitively computationally expensive, as the time step scales as 

1 /2 

5t (X kfi for collisions to be modeled effectively. Previous simulations have shown that increasing k^i does 
not significantly alter physical results [15j. We use a time step of St = 1.0 x 10~^r and damping coefficients 
7n = 50r~^ and = O-O- All measurements are expressed in terms of d and r. We make use of both 
cartesian coordinates (x, z) and cylindrical coordinates (r, 0, z) where = + ?/^; gravity is set to ldT~^ 
and points in the negative z direction. 

The simulations were carried out on MIT's Applied Mathematics Computational Laboratory, a Beowulf 
cluster consisting of 16 nodes each with a dual processor. During the simulations, snapshots of all particle 
positions were saved every 2r, corresponding to 20000 integration timesteps. The LAMMPS code is written 
in C++ and can be run on any number of processors, by decomposing the computational domain into 
cuboidal subdomains of equal size. Interactions between particles in neighboring domains are handled using 
message passing. For problems where the particles are split evenly between the subdomains, the LAMMPS 
code scales very well, and doubling the number of processors can frequently result in almost a doubling of 
speed. However, the Couette geometry causes some subdomains to have more particles than others, leading 
to poor load-balancing. Most of the simulations were therefore run on four processors, with the domains 
split by the planes x = and = to achieve optimal load-balancing; these simulations took approximately 
four weeks to complete. 

4.2 The silo geometry 

We considered a quasi-2D silo with plane walls at x = ±75(i, z = with friction coefficient ji^ = 0.5, and 
made the y direction periodic with width Sd. To generate an initial packing, 90, 000 spherical particles with 
contact friction coefficient jHc = 0.5 were poured in from a height of 2: = 130d at a rate of 378r~^. After all 
particles are poured in at t = 238r, the simulation is run for an additional 112r in order for the particles to 
settle. After this process has taken place, the particles in the silo come to a height of approximately z = 62.2d. 
To initiate drainage, an orifice in the base of the container is opened up over the range —3d<x<3d, and 
the particles are allowed to fall out under gravity; Figure [3l^a) shows a typical simulation snapshot during 
drainage. 
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As noted in the previous section, snapshots of the particle positions are saved every 2r. We cohected 282 
snapshots, and made use of this information to construct velocity cross sections. A particle with coordinates 

at the nth timestep and x^+i at the (n + l)th timestep makes a velocity contribution of (x^+i — x^)/At 
at the point with coordinates (x^ +Xn+i)/2. This data can then be appropriately binned to create a velocity 
profile; we considered bins of size d in the x direction, and created velocity profiles for different vertical slices 
\z - Zs\ < d/2. 

Since the SFR makes predictions about the velocity profile during steady flow, we choose a time interval 
ti < t < t2 over which the velocity field is approximately constant. Choosing this interval requires some 
care, since if ti is too small then initial transients in the velocity profile can have an effect, and if t2 is too 
large, then the free surface can have an influence. For the results reported here, we chose ti = 120r and 
t2 = 200r. 

Figure EJb) compares the SFR predictions for this environment to the DEM simulation. The displayed 
simulation data uses a particle contact friction of /ic = 0.3. Since the typical range of Ls from velocity corre- 
lations in simulations [30j and experiments [14J is 3d to 5(i, we choose Ls = M to generate the approximate 
SFR solution. We emphasize that this parameter is not fitted. In this geometry, the slip-lines are symmetric 
about the drift direction causing both Di and D2 to become independent of the internal friction. In prior 
simulations we have found that that particle contact friction has some effect on the flow [31] , and analogously 
the internal friction should play some role in the determination of b. Here, the loss of friction dependence 
comes from our simplification that D2 is isotropic. A less simple but more precise definition for D2 would 
anisotropically skew the spot diffusion tensor as a function of internal friction: the slip-lines, which we model 
as the directions along which a spot can move (roughly), intersect at a wider angle as internal friction is 
increased. 

Even so, our simple model captures many of the features of the flow and accurately portrays the dominant 
behavior. The downward velocity, especially at z = lOd, strongly matches the predicted Gaussian. Perhaps 
a more global demonstration of the underlying stochastic behavior in the SFR is evident in Figure [31(c), 
where a linear relationship can be seen between the mean square width of Vz and the height, indicating that 
the system variables are undergoing a type of diffusive scaling. The SFR solution also predicts this linear 
relationship and, in particular, that the slope should equal 2b = Lg. The agreement shown in Figure 7(b) 
for such a typical Ls value is a strong indicator that the role of the correlation length in the flow has been 
properly accounted for in the SFR. 

The diffusive velocity profile near the orifice of a wide silo has also been observed in a number of ex- 
periments ^25l [371 O [32] and DEM simulations [301 ISl and was first predicted qualitatively by kinematic 
models [TTl [24l [27]. Such models, however, lack a plausible microscopic basis [4 or any connection with 
mechanics and cannot be applied to any other geometry, other than the wide silo. 

4.3 The Couette geometry 

For the Couette geometry, we considered five different interparticle friction coefficients, /ic = 0.1, 0.3, 0.5, 0.7, 0.9, 
and for each value an initial packing was generated using a process similar to that for the silo. We consider 
a large cylindrical container with a side wall at r = 64(i with friction coefficient /j.^ = /^c, and a base at 
z = with friction coefficient jiiw = 0. For each simulation, 160000 particles are poured in from a height 
of z = 4:Sd at a rate of 4848r~^. After all particles are introduced at t = 33r, the simulation is run for an 
additional time period of 322r to allow the particles to settle. After this process has taken place, the initial 
packings are approximately 11.5d thick. Packings with /ic = 0.9 are approximately 2% thicker than those 
with /ic = 0.1. 

The shearing simulations take place in an annulus between two radii, rin and rout- A rough outer wall 
is created by freezing all particles which have r > rout; any forces or torques that these particles experience 
are zeroed during each integration step. Similarly, all particles between rin — M and rin are forced to rotate 
with a constant angular velocity uo around r = 0. Particles with r < rin — 4(i are deleted from the simulation 
and an extra cylindrical wall with friction coefficient /i^^ = /ic is introduced at r = rin — 4(i to prevent stray 
surface particles from falling out of the shearing region. A typical run is shown in Figure [It^ a). 
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(b) (c) 

Figure 3: (a) A typical snapshot of the silo system during drainage, taken at t = 60r. The colored bands 
are initially spaced lOd apart, and highlight the deformations that occur during flow (b) Mean downward 
velocity profiles at two different heights in a wide silo. DEM results (solid lines) are plotted against the SFR 
predictions (dashed lines) at heights of z = lOd, 30d; (c) Mean square width of the downward velocity across 
horizontal cross-sections in the DEM (points), compared to the SFR (dashed line) 



From the five initial packings, we carried out eight different shear cell simulations to investigate the effects 
of friction, angular velocity, and inner wall radius. To investigate the effect of friction, five runs were carried 
out with uj = O.Olr"^ and rin = 40(i, for /ic = 0.1,0.3,0.5,0.7,0.9. To examine the effect of the inner wall 
radius, an additional run with rin = 30d and rout = SOd was carried out, with /ic = 0.3 and uo = O.Olr"^ kept 
constant. To look at the effect of angular velocity, a further two runs with co = 0.05r~^ and co = 0.2r~^ were 
carried out, with jHc = 0.5 and rin = 40(i kept constant. For each simulation, we collected 561 snapshots. For 
the runs where rin = 40(i, approximately 108350 particles were simulated, corresponding to 2.0Gb of data. 
For the run with rin = 30(i, 88657 particles were simulated, corresponding to 1.7Gb of data. 

The simulation results show a good empirical agreement with previous experimental work on shear cells 
[IHl Ell Uni El [23]. In all cases, we see an angular velocity profile that falls off exponentially from the inner 
cylinder, with a half-width on the order of several particle diameters. Near the inner wall, the flow deviates 
from exponential, which is an effect seen in some prior studies but is more dramatic here. Following methods 
similar to those used in silo simulation, we used the snapshots to construct an angular velocity profile. We 
used bins of size d/2 in the radial direction, and we looked at velocity profiles in different vertical slices 

< Z < ^high- 

Since the simulation geometry is rotationally symmetric, our angular velocity profile can most generally 
be a function of r, z and t. Ideally, we hope that co is primarily a function of r, with only a very weak 
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Figure 4: (a) A typical snapshot of the system during shearing based on the simulation with rin = 40(i, 
CO = O.Olr"^, and /ic = 0.5. Dark blue particles (with r > rout) are frozen during the simulation, while dark 
red particles (with r < rin) are rotated with angular velocity co. The particles between rin and rout undergo 
shearing. The colored bands of grey, cream, and cyan were initially radial, and in this snapshot, after fifty 
frames, the deformations can be clearly seen, (b) SFR solutions using wall hypothesis 1 (dashed line) and 
wall hypothesis 2 (dot-dashed line) in the annular Couette geometry (rin = 40(i, rout = 60(i) are compared 
to the simulation results (solid line). Both use the same parameters that were used in the silo geometry 
{Ls = ^d^ lie = 0.3) and the bead internal friction angle is set at the typical value 25^. 



dependence on z and t, but we began by studying the effects of these other variables. To determine the 
dependence on time, the velocity profiles were calculated over many different time intervals. As would be 
expected, the simulation had to be run for small amount of time before the velocity profile would form; this 
happens on a time scale of roughly 50r, and the results suggest a longer time is needed for the cases with low 
friction. However, the data also shows time-dependent effect happening on a longer scale: as the shearing 
takes place, there is a small but consistent migration of particles away from the rotating wall, which has a 
small effect on the velocity profile. This effect does eventually appear to saturate, but because of this, we 
chose to discard the simulation data for t < 500r and calculate velocity profiles based on the time window 
500r < t < llOOr. 

To investigate the angular velocity dependence on the height, we calculated the velocity profiles in five 
different slices Zh < z < Zh -\- d for Zh = d^ Sd, Sd, Td, 9d. Near the inner rotating wall, the velocity profiles 
show very little dependence on height. However, in the slow-moving region close to the fixed outer wall, 
large differences can be seen, with particles in the lowest slice moving approximately 30% slower than those 
in the central slice, and those in the top slice moving approximate 30% faster. The three central slices show 
differences of at most 10%, and we therefore chose to use the range 3d < z < Sd. 

The SFR treats the correlation length as a material property independent of the flow geometry or other 
state variables. To see how well this notion is upheld, we solve the SFR in the annular Couette geometry 
using the same correlation length {Lg = M) that was used in the displayed silo prediction. Figure [3l It is 
then compared to a simulation of annular flow which uses the same grain properties (/ic = 0.3). 

Results from Figure [4](b) show that regardless of the wall hypothesis, the SFR prediction captures the 
same qualitative features of the simulation. The SFR and simulation both predict a flatter range near the 
inner wall, followed by exponential decay. Near the inner wall, it does appear that wall hypothesis 1 ("no 
slip condition" for spots) gives a closer match to the simulation. 

The SFR, when applied to the annular flow geometry, does predict a slight dependence of the flow on the 
internal friction. We emphasize that internal friction is not the same quantity as particle contact friction 
/ic5 though we believe if the contact friction is increased, inevitably, the internal friction must be as well. 
Spherical grains almost always have internal friction angles in the range (f) = 20° to ^ = 30°, so to represent 
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Figure 5: Velocity profiles for five different values of /ic, where rin = 40(i and uo = O.Olr 
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Table 1: (Left) Simulation: Half- widths of the shearing velocity profiles for different values of /ic, calculated 
by fitting the functional form /(r) = a — (log2)r/6 to Xogv/v^ over the range Abd < r < bSd. (Right) SFR: 
Fits the predictions to the same form and uses Lg = 3d. 



this range as best as possible, we simulated flows varying the particle contact friction from jHc = 0.1 to 
/ic = 0.9. Figure [5] displays velocity profiles for five different values of friction. In the semi-log format, 
profiles appear almost linear over the range 45(i < r < 5Sd indicating a very good fit to an exponential 
model of velocity. 

Since the curves in the figure are very close, and exhibit some experimental noise, it is difficult to discern 
any small differences in the widths of the velocity profiles. However, table [T] shows the results of applying 
linear regression to extract a half- width for each velocity profile. We see differences on the order of 5%, 
roughly in line with the SFR. More importantly, the trend of increasing flow width with increasing friction 
is seen in both. 

When the inner wall radius is decreased. Figure [6] indicates that the shear band shrinks but the decay 
behavior in the tail changes only minimally. The SFR predicts this qualitative trend as well, but significantly 
underestimates the size of the shear-band decrease. 

In agreement with past work on Couette flow [TSl [5] , we too find that the normalized flow profile is roughly 
unaffected by the shearing rate (see Figure [7j). As previously discussed, this behavior is in agreement with 
the SFR, which always permits flow fields to be multiplied by a constant. 



5 Conclusion 

We have presented the crucial principles motivating the Stochastic Flow Rule and have assessed its validity by 
checking analytical predictions for silo and annular Couette flow against discrete-element simulations. Using 
the same parameters for both cases, without any fitting, the SFR manages good predictions for these two 
very different flow geometries, which it seems cannot be described, even qualitatively, by any other model. 
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Figure 6: DEM simulation results (solid lines) for the flow profile when the inner wall radius is 30d (blue) 
and 40d (green) compared to SFR predictions (dashed lines) which use the same color scheme. For solidarity 
with the silo flow predictions, the SFR solutions use Ls = M. 




Figure 7: Velocity profiles for three different angular velocities, with /ic = 0.5 and rin = 40(i. The time 
windows over which the velocities are computed are scaled according to the angular velocity. 
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The model was also shown to capture the "diffusive" type flow properties unique to granular materials such as 
Gaussian downward velocity in the draining silo and exponentially decaying velocity in the annular Couette 
ceh. 

Our simulations also indicate that the slight changes in flow brought on by varying the inter-particle 
contact friction in the annular flow geometry match the trends the SFR predicts when the internal friction 
angle is approrpiately varied. The trend is also captured when the inner wall radius is varied, though the 
quantitative agreement is not as strong. In agreement with past studies on annular Couette flow, and in 
validation of one of the flrst principles behind the SFR, we flnd in our simulations that the flow rate does 
not signiflcantly affect the normalized flow proflle over a signiflcant range of rates. 

These results, taken together with our prior work on spot-driven particle simulations [30], suggest that 
the SFR provides a general paradigm for multiscale modeling of granular flow or, more generally perhaps, 
plastic deformation of other amorphous materials. Whenever the stress state is near yielding our physical 
picture of flow is that spots of fluidization perform random walks along slip-lines biased by local stress 
imbalances. With some exceptions as outlined in [T4] (so-called "slip-line admissible" geometries) we believe 
the SFR is a highly general law for 2D dense granular flow. 
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